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ABSTRACT 

Context. There is a need to improve the fidelity of SPH simulations of self-gravitating gas dynamics. 

Aims. We remind users of SPH that, if smoothing lengths are adjusted so as to keep the number of neighbours, N, in the range 
± AN NmB , the tolerance, AJV NEIB , should be set to zero, as first noted by Nelson & Papaloizou. We point out that this is a very 
straightforward and computationally inexpensive constraint to implement. 

Methods. We demonstrate this by simulating acoustic oscillations of a self-gravitating isentropic monatomic gas-sphere (cf. Lucy), 
using N Tm ~ 6, 000 particles and N NUB = 50. 

Results. We show that there is a marked reduction in the rates of numerical dissipation and diffusion as AJV NEIB is reduced from 10 to 
zero. Moreover this reduction incurs a very small computational overhead. 

Conclusions. We propose that this should become a standard test for codes used in simulating star formation. It is a highly relevant 
test, because pressure waves generated by the switch from approximate isothermality to approximate adiabaticity play a critical role 
in the fragmentation of collapsing prestellar cores. Since many SPH simulations in the literature use A/^g = 50 and AN NE1B > 10, 
their results must be viewed with caution. 
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1. Introduction 

Numerical simulations play an increasingly important role in the 
study of star formation, and of other non-linear phenomena in- 
volving self-gravitating gas dynamics, for example the develop- 
ment of cosmological structure, galaxy formation and stellar col- 
lisions. Since these problems involve complex three-dimensional 
configurations and/or large ranges of density, Smoothed Particle 
Hydrodynamics (SPH; e.g. Monaghan 1992) is well suited to 
such simulations. The principal alternative is to use an Adaptive 
Mesh Refinement (AMR) Finite Difference code (e.g. Truelove 
et al. 1998), although this is in general more expensive computa- 
tionally. However, ultimately neither method is useful, unless it 
can be shown that the results are converged, and are not compro- 
mised by numerical artefact. This requires test problems with 
known analytic, or semi-analytic, solutions, which involve the 
same basic physical phenomena as occur during star formation. 

We propose such a test here, namely acoustic oscillations of a 
self-gravitating, isentropic, monatomic gas-sphere in the funda- 
mental radial mode. This test was originally performed by Lucy 
(1977) in his seminal paper introducing SPH, and therefore we 
shall refer to it as the Lucy test. It has been performed subse- 
quently (e.g. Steinmetz & Midler 1993; Nelson & Papaloizou 
1994), but infrequently. It is an appropriate test because it mea- 
sures (i) the level of dissipation associated with artificial viscos- 
ity, in the absence of shocks; (ii) the extent to which transients, 
due to the discrete nature of particles (or cells), remove energy 
from genuine modes and transfer it to other spurious modes (i.e. 
numerical diffusion); (iii) the ability of the code to model acous- 
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tic oscillations, and in particular adiabatic bounces; and (iv) the 
ability of the code to deal with free (or nearly free) boundaries. 

Point (iii) is particularly important because it seems that 
collapsing prestellar cores are most prone to fragment at the 
stage when the gas switches from being approximately isother- 
mal to being approximately adiabatic (e.g. Boss et al. 2000). 
Fragmentation at this juncture is presumably due to interactions 
between the complex system of pressure waves which is gener- 
ated by adiabatic bounces in a converging but disordered inflow, 
and it is therefore essential that spurious waves are not being 
generated. In this context it is worth noting that the isentropic 
assumption is not strictly the same as adiabaticity, and is made 
here for analytic convenience rather than realism. In simulations 
where shocks develop, artificial viscosity must be incorporated, 
in which case, either the resulting energy dissipation must be in- 
cluded in the energy equation, or - if the thermal timescale is 
sufficiently short, and the main concern is not with the detailed 
structure of the shock front - a barotropic equation of state may 
be invoked. There are no shocks in the present simulation, and 
therefore the rate of energy dissipation due to artificial viscosity 
is low, but not negligible. The isentropic assumption then simply 
implies that the small amount of energy dissipated by artificial 
viscosity is spirited away by some unspecificed process. 

It is also appropriate to point out that, although the gas in 
prestellar cores is largely molecular, it behaves as a monatomic 
gas (i.e. effective adiabatic exponent y = 5/3) until the temper- 
ature rises above ~ 100 K. At lower temperatures the rotational 
degrees of freedom are not significantly excited, since para-H, 
has its first excited level (J = 2) at £ B (512K), and ortho-H, 
has its first excited level (J = 3) at fc B (854K) (e.g. Black & 
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Bodenheimer 1975). Neglect of this fact can lead to artificially 
enhanced fragmentation of collapsing prestellar cores, since with 
y = 5/3 the Jeans mass increases quite rapidly with increasing 
density in the adiabatic regime (M JEANS oc p 1 / 2 ), whereas with 
y = 7/5 it increases much more slowly (M JEANS oc p'^ 10 ). 

In Section [2] we describe the numerical code used and how 
the test is initiated. In Section [3] we present and discuss the re- 
sults. In Section|4]we summarise our conclusions. 



Table 1. Decay of the fundamental mode. The first column gives 
AAf NEIB . The second column gives the e-folding time for the am- 
plitude of the fundamental mode, T A , as a multiple of its period, 
P . The third column gives the net oscillation energy left after 10 
periods, £ 10 , as a fraction of the initial oscillation energy, £ . The 
fourth column gives the simulation time to evolve the oscillating 
gas-sphere for 10 periods on eight 2.2Ghz Opteron CPUs each 
with 8GB memory, f 10 . 



2. Numerical details 

The code we have used - dragon (Goodwin et al. 2004) - is 
a standard SPH code, using a second-order Runge-Kutta inte- 
gration scheme, multiple particle timesteps, adaptive smoothing 
lengths, time-dependent artificial viscosity (Morris & Monaghan 
1997), and an octal spatial-decomposition tree (Barnes & Hut 
1986) for calculating gravitational accelerations and construct- 
ing neighbour lists. In calculating gravitational accelerations, 
tree-cells are opened if they present an angle greater than 
9 CRn = 0.5, and quadrupole moments are included. The smooth- 
ing length of a particle is adjusted so that its neighbour list con- 
tains N = Af NEIB ± &N mrB other particlefl 

As an illustration of the Lucy test, we demonstrate how much 
the performance of this code improves as &N NElR is reduced, and 
how marginal the extra computational cost is. The equilibrium 
isentropic sphere is set up by placing equal-mass particles on an 
hexagonal close-packed array within a sphere of radius unity; 
then stretching this uniform-density sphere radially to reproduce 
the density profile of a poly trope with exponent 5/3 (or equiva- 
lently index 3/2); and finally relaxing the system by evolving the 
particle positions using the SPH code, until the net kinetic energy 
falls to a very small value (relative to the magnitude of the gravi- 
tational potential energy). In the results presented here the sphere 
is represented by A/ T0T = 5, 895 particles, and we set N NEm = 50. 
The parameters in the time-dependent artificial viscosity (Morris 
& Monaghan 1997) are a t = 0.1, C, = 0.2,/? = 2a. To excite 
the fundamental mode, each particle is displaced radially from 
its equilibrium radius r to a new radius r' = r [1 + A^ir/R)], and 
then released from rest. Here <;(s) (0 < s < 1) is the eigenfunc- 
tion of the fundamental radial modfl A is the initial amplitude 
of the oscillation, which we set to A = 0.1. 

3. Results and Discussion 

Figures 1 and 2 show oscillations simulated with AN NElB = 
0, 2, 5, and 10. In Fig. 1, the quantity plotted is the mean x- 

1 Some SPH codes (e.g. Price & Monaghan 2004) do not specify the 
number of neighbours, but instead, for each particle, iterate around a 
loop, 



until fractional changes in p ; drop below a user-defined tolerance, e. 
Here h a is a constant of order unity, and the summation is over all neigh- 
bours j for which - r.| < 2 h j (i.e. all particles within the smoothing 
kernel of particle (')■ Provided e is sufficiently small, this is statistically 
equivalent to setting Af NEIB = 32nh 3 g /3 and AN NEm = 0. 

2 The eigenfunction, £(s), for the fundamental radial mode of a self- 
gravitating isentropic monatomic gas-sphere was very kindly supplied, 
in the form of a dense look-up table, by Alfred Gautschy. 
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(the mean being taken over all the particles, with equal weight- 
ing), as a function of time (normalised to the central freefall time, 
t w , in the unperturbed equilbrium state). In Fig. 2, the quantity 
plotted is the total kinetic energy, normalised to the magnitude 
of the gravitational potential energy in the equilibrium state (i.e. 
< K/\S1 \, where 7C is the total kinetic energy, Q. is the total grav- 
itational potential energy, and subscript refers to values in the 
unperturbed equilibrium state), again as a function of time (nor- 
malised to fpp). Some decay parameters are recorded in Table 1. 

In addition to the fundamental mode, some additional modes 
are unintentionally excited from the outset. This is because, 
following relaxation, the equilibrium state of an isentropic 
monatomic sphere is not modelled exactly, due to the smooth- 
ing inherent in SPH. In particular, the density is underestimated 
near the centre and near the boundary. (This can be improved by 
increasing N NEm and N TOT /N NEm , but that requires extra compu- 
tational resource.) Furthermore, the fundamental mode is excited 
with finite amplitude, but the eigenfunction is derived on the as- 
sumption of infinitesimal amplitude. (This can be moderated by 
adopting a lower value for A, but the test is not then relevant to 
real simulations of star formation, where ultimately we are con- 
cerned with non-linear perturbations.) 

Setting aside the effect of other modes present in the initial 
conditions, the subsequent decay of the fundamental mode is in 
general due to two effects. First, the oscillation energy may be 
dissipated by artificial viscosity. The dissipation of energy due to 
artificial viscosity occurs wherever two neighbouring SPH par- 
ticles approach one another. Second, the oscillation energy may 
be transferred to other modes. This occurs wherever local den- 
sity or velocity fluctuations are created by the discrete nature of 
the SPH particles, or by the 'peculiar velocities' of individual 
SPH particles. Both effects occur more rapidly for larger values 
of AN NEm . (They also occur more rapidly for smaller values of 
N Tar and smaller values of Af NEIB , but these parameters are not 
varied here.) 

When AAf NEIB = 0, the neighbour list of an SPH parti- 
cle changes very seldom, and - when it does - very little. 
Therefore the acceleration experienced by the particle varies 
very slowly, and the velocity field is very smooth. The upshot 
is that neighbouring particles only approach one another very 
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Fig. 1. The mean x-displacement, x (where the mean is over all particles, with equal weighting), against time, t (in units of the 
freefall time in the equilibrium state, t w ). The results are displaced vertically by Ax to fit them all on one plot. Reading from the top, 
(a) AN NE[B = 0, Ax = 0; (b) AN NEm = 2, Ax = -0.05; (c) AN NE[B = 5, Ax = -0.10; (d) AAf NEIB = 10, Ax = -0.15 . 



slowly, and the rate of dissipation due to artificial viscosity is 
low. Notwithstanding the slow rate of dissipation, there are small 
fluctuations in density and velocity, and these feed energy into 
other modes, so that the fundamental mode decays (see Figs. 1 
& 2). 

As AN NEiB is increased, the neighbour list of an SPH particle 
changes more frequently, and more abruptly. Therefore the ac- 
celeration experienced by the particle varies in a more idiosyn- 
cratic manner, and the velocity and density fields are more noisy. 
The upshot is that neighbouring particles often approach one 
another more rapidly, and the rate of dissipation due to artifi- 
cial viscosity is therefore higher. In addition, the noisy velocity 
and density fields are very effective at feeding energy into other 
modes, so that the fundamental mode decays more rapidly (see 
Figs. 1 & 2). 



2AJV NEIB 



In principle, the number of neighbours can change by 
at each timestep, from N NEm - AN KElB to N NEm + AN NElB 
or vice versa. Thus with A/ NEIB = 50 and AAf NEIB = 10, the number 
of neighbours can change from 40 to 60 or vice versa. In prac- 
tice such large changes are unlikely, but it is still the case that 
increasing AN NEm results in increased fluctuations in the neigh- 
bour lists. In particular, particles in condensing regions tend to 



whilst particles in expanding regions tend to have N fluctuating 
between ~ A/ NEIB and ~ (Af NEIB - AAf NEIB ). 

In Table 1, we record, for each value of AN mm (column 1), 
the e-folding time of the amplitude of the fundamental mode, 
T A , in terms of its period, P = 3.7f FF (column 2); the oscillation 
energy left after ten periods, £ 10 , as a fraction of its initial value, 
£ (column 3); and the simulation time to evolve the oscillating 
gas-sphere for 10 periods on eight 2.2Ghz Opteron CPUs each 
with 8GB memory (column 4). The oscillation energy is given 
by 



•7C + (T-T ) + (Q-n ), 



(2) 



have N fluctuating between 



N NE[B and 



where T is the thermal energy, and again the subscript indi- 
cates the unperturbed equilibrium state. 

We note that the oscillation energy decays rather slowly due 
to dissipation, even with AN^ = 10. This is because the oscil- 
lations have low amplitude, and therefore the relative velocities 
between neighbouring particles are always very subsonic. Not 
only are the initial amplitudes low, but in the cases with high 
AAffjHu the amplitudes decay rapidly due to numerical diffusion. 
In other words, when AN^ EXB is high, the rate of dissipation is 
reduced because diffusion rapidly spreads the oscillation energy 
amongst many different modes and thereby reduces even further 
the relative velocities between neighbouring particles. This is re- 
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Fig. 2. The total kinetic energy, 'K (normalised to the magitude of the self-gravitational potential energy in the equilibrium state, 
|O 1), against time, t (in units of the freefall time in the equilibrium state, f FF ). The results are displaced vertically by A'K to fit them 
all on one plot. Reading from the top, (a) AJV NEIB = 0, A'K = 0; (b) AJV NEIB = 2, A'K = -0.003; (c) A7V NEIB = 5, A'K = -0.006; 
(d) AJV NEIB = 10, A'K = -0.009 . 

liable than those obtained with AN NEIB = 10, and significantly 
more reliable even than those obtained with AN NEm = 2, both 
in terms of having less dissipation and - more importantly - in 
terms of having less numerical diffusion. 

Setting AN NEIB = also requires little extra computation. For 
example, to follow 10 oscillation periods with AN NE[B = takes 
only 17% longer than with A,/V NEIB = 10. Moreover, 7% of this 
increase is due to the fact that with AN NElB = the sphere contin- 
ues to oscillate with a significant amplitude after 10 periods, and 
therefore the timestep remains short. When allowance is made 
for this, the real cost of reducing AN NElB from 10 to is only a 
10% increase in computation. 

Therefore our principal conclusions are (i) that AN NE , B = 
should be the default option for SPH codes which adapt smooth- 
ing lengths in this way; and (ii) that the Lucy test provides a very 
useful way of evaluating the fidelity of codes used in simulations 
of star formation. 

We emphasise that we have not set out to reproduce as accu- 
rately or economically as possible acoustic oscillations of a self- 
gravitating isentropic sphere in the fundamental mode. We have 
simply demonstrated how the results produced using a standard 
SPH code, with a modest number of particles (N TQT = 5,895) 
depend on AN NEiB . There are adjustments to SPH which will im- 



flected in the results presented in Fig. 2 and the third column of 
Table 1 . Because the decay of the fundamental mode is largely 
due to numerical diffusion, the oscillation energy only falls by 
a few percent after ten periods (see the third column of Table 
1), whereas the amplitude of the variation in kinetic energy falls 
much more rapidly. The variation in kinetic energy eventually 
disappears, not because the kinetic energy itself disappears, but 
because numerical diffusion transfers oscillation energy to other 
modes with different periods and different phases. Consequently 
the oscillation energy becomes thermalised, and 7C is finite but 
approximately constant. 

The simulations presented here have been continued for 100 
periods. In this limit, there are so many modes excited, with so 
many different phases, that the gas-spheres become virialised 
with 

27C + 2T + Q ; (3) 
% is still finite. 

4. Conclusions 

From the plots in Figs. 1 and 2, and the above discussion, we 
infer that the results obtained with AN n „ B = are far more re- 
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prove (extend) the timescales for dissipation and numerical dif- 
fusion in the present simulation. For example, using standard ar- 
tificial viscosity with a = and j3 — 0.1 increases the e-folding 
time for the fundamental mode to ~ 60 oscillation periods, but 
at the same time compromises the shock-capturing ability of the 
code so that it can not then be used for simulations in which 
shocks are likely. Similarly, the e-folding time for the fundamen- 
tal mode can be extended by increasing N NmB or M TOT /N , but 
this must inevitably be at the expense of resources which are 
needed elsewhere, viz. to maximise the extent and/or duration of 
a simulation. 

The useful duration of the present simulation can be iden- 
tified with the e-folding time of the fundamental mode, which 
withA,/V NEIB = is 13.6 oscillation periods, but with AAf NEIB = 10 
is only 3.1 oscillation periodfl 

Adaptive Mesh Refinement codes used for simulating star 
formation should also aim to reproduce or improve upon the re- 
sults produced here, using comparable computational resources. 
In addition, they should do so for a gas-sphere which moves at 
constant velocity relative to the underlying Cartesian grid, in or- 
der to match the Lagrangian advantages of SPH. 
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